
/**
 输出:
 1. UTC改为BJ
 2. 加星时
 2024-09-21T10:00:00 <==> 54381600
 3. 方位角阈值 ±10°以外
 4. 高度角阈值: 30-80°
 */

#include <stdio.h>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <string>
#include <algorithm>
#include "ADefine.h"
#include "ATimeSpace.h"
#include "OrbitDef.h"
#include "SimStar.h"
#include "ExtractWCS.h"

using namespace boost;
using namespace boost::posix_time;
using namespace boost::filesystem;
using namespace AstroUtil;

void Bjt2DateTime(int st, int& date, int& time) {
    int st0 = 54381600;
    ptime bjt0 = from_iso_extended_string("2024-09-21T18:00:00");
    ptime bjt = bjt0 + seconds(st - st0);
    ptime::date_type day = bjt.date();
    ptime::time_duration_type hms = bjt.time_of_day();
    int YY = day.year();
    int MM = day.month().as_number();
    int DD = day.day();
    int hh = hms.hours();
    int mm = hms.minutes();
    int ss = hms.seconds();
    date = (YY * 100 + MM) * 100 + DD;
    time = (hh * 100 + mm) * 100 + ss;
}

double Bjt2MjdUTC(int st) {
    int st0 = 54381600;
    ptime utc0 = from_iso_extended_string("2024-09-21T10:00:00");
    ptime utc = utc0 + seconds(st - st0);

    return utc.date().modjulian_day() + utc.time_of_day().total_milliseconds() / 86400000.0;
}

void VerifyPoint() {
    ATimeSpace ats;
    ats.SetSite(116.26663, 40.06856, 88, -8);

    int st  = 57156896;
    int st0 = 54381600;
    ptime bjt0 = from_iso_extended_string("2024-09-21T18:00:00");
    ptime bjt = bjt0 + seconds(st - st0);
    ptime utc = bjt - hours(8);
    ptime::date_type day = utc.date();
    ptime::time_duration_type hms = utc.time_of_day();
    int YY = day.year();
    int MM = day.month().as_number();
    int DD = day.day();
    ats.SetUTC(YY, MM, DD, (hms.total_seconds() + 0.276) / 86400.0);
    printf ("%s\n", to_iso_extended_string(utc).c_str());
    printf ("epoch: %lf\n", ats.Epoch());
    double ra2k(342.773688), dec2k(-11.052835);
    double raapp, decapp;

    ats.EqTransfer(ra2k * AU_D2R, dec2k * AU_D2R, raapp, decapp);
    printf ("%lf  %lf\n", raapp * AU_R2D, decapp * AU_R2D);
}

bool ClipOrbit(int argc, char **argv) {
    if (argc == 1) return false;

    OrbPtrVec orbs; // 轨迹集合
    OrbitPtr orbOne;  // 单轨迹

    FILE *fp = fopen(argv[1], "r");
    char line[1024];

    int code, codeOld(-1), st, stOld(-100);
    double dur, vel, azi, alt, ra, dec;
    double ra2k, dec2k, mjd;
    char bjtstr[40];

    ATimeSpace ats;
    ats.SetSite(116.26663, 40.06856, 88, -8);

    while (fgets(line, 1024, fp)) {
        // 解析
        sscanf(line, "%d %s %d %lf %lf %lf %lf %lf %lf",
            &code, bjtstr, &st, &dur, &azi, &alt, &vel, &ra, &dec);

        if (alt < 30.0 || alt > 80.0) continue;
        if (azi < 10.0 || azi > 350.0) continue;

        if (code != codeOld || st < stOld || (st - stOld) > 60) {
            codeOld = code;

            if (orbOne.unique()) {
                if (orbOne->Size() >= 10) orbs.push_back(orbOne);
                else printf ("%6d is abandoned\n", orbOne->code);
            }

            orbOne.reset(new Orbit);
            orbOne->code = code;
            Bjt2DateTime(st, orbOne->date, orbOne->time);
        }
        stOld = st;

        // 历元转换: 当前-->J2000
        ats.SetMJD(Bjt2MjdUTC(st));
        ats.EqReTransfer(ra * AU_D2R, dec * AU_D2R, ra2k, dec2k);
        // 历元转换: 结束

        OrbitNode node;
        strcpy(node.bjtstr, bjtstr);
        node.tmSat = st;
        node.azi = azi;
        node.alt = alt;
        // node.ra  = ra;
        // node.dec = dec;
        node.ra  = ra2k * AU_R2D;
        node.dec = dec2k * AU_R2D;
        orbOne->nodes.push_back(node);
    }
    fclose(fp);
    fp = NULL;

    // 排序
    std::stable_sort(orbs.begin(), orbs.end(), [](const OrbitPtr& o1, const OrbitPtr& o2) {
        return o1->nodes[0].tmSat <= o2->nodes[0].tmSat;
    });
    // 删除重复弧段
    int time = -1;
    code = -1;
    for (auto&it : orbs) {
        if (code == it->code && time == it->time) it->repeated = true;
        code = it->code;
        time = it->time;
    }

    // 输出总文件
    printf ("output >> total >> begin\n");

    int date = 0;
    for (auto &it : orbs) {
        if (it->repeated) continue;

        if (it->date != date) {
            if (fp) {
                fclose(fp);
                fp = NULL;
            }
            date = it->date;
            printf ("resolving : %d\n", date);
        }
        if (!fp) {
            // 创建目录
            path pathName((boost::format("total_%d.txt") % date).str());
            // 打开文件
            fp = fopen(pathName.c_str(), "w");
        }
        printf ("%6d  %d\n", it->code, it->nodes[0].tmSat);
        fprintf (fp, "%s\n", std::string(100, '-').c_str());
        it->Output(fp);
    }
    if (fp) fclose(fp);
    printf ("output >> total >> over\n");
    // 输出单轨文件
    printf ("output >> single orbit >> begin\n");
    path pathDir;
    date = 0;
    for (auto &it : orbs) {
        if (it->repeated) continue;

        if (it->date != date) {
            date = it->date;
            pathDir = path((boost::format("WY%d") % date).str());
            if (!exists(pathDir)) create_directories(pathDir);
        }
        boost::format fmt("%dT%06d_%d.txt");
        fmt % it->date % it->time % it->code;
        path pathName(pathDir);
        pathName /= fmt.str();
        printf ("\t %s\n", fmt.str().c_str());

        fp = fopen(pathName.c_str(), "w");
        it->Output(fp);
        fclose(fp);
    }
    printf ("output >> single orbit >> over\n");

    printf ("Over\n");

    return true;
}

int main(int argc, char **argv) {
    // GenerateSimStar("2024-10-25T00:00:00", "2024-10-25T05:00:00", 5);

    ClipOrbit(argc, argv);

    // if (argc < 2) return -1;
    // ExtractWCS(argv[1]);

    return 0;
}
